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ABSTRACT 


As linear programming techniques find applications in more 
diverse fields, the problem of solution time becomes increasingly 
important. A variation of the revised simplex algorithm, in which 
the constraints are added in a step-by-step fashion, is investigated 
as a potentially faster solution technique. A computational pro- 
cedure, coded for the IBM 360 computer, is developed to compare this 
algorithm with the standard two-phase revised simplex algorithm. 

A limited number of problems, including several randomly generated 
problems, is solved by each of the two methods. The resulting com- 
parison of solution times indicates that a significant improvement 
is obtained by the use of the procedure of step-by-step addition of 


constraints. 
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ls Introduction 

Increased utilization of linear programming techniques has led to 
the formulation of problems of sufficient size to tax the storage 
capabilities of many computers now in use. The problem of storage 
capacity may be alleviated by auxiliary storage, with the accompanied 
requirement of access delays which increase solution time. Even with 
sufficient storage and a coded procedure large enough to accommodate 
a problem, the computation time may be such that solution costs ap- 
proach a budget limitation. For example, the CEIR LP/90 program for 
the IBM 7090 computer, which can accommodate 512 constraints and an 
unlimited number of variables, or the Philco Corporation LP 2000 
System for the S=2000 computer, which can accommodate 2500 constraints 
and an unlimited number of variables, might take several hours to 
solve a large problem at a cost of several thousand dollars. 

Improvement in the solution time might be accomplished by im- 
provement of the hardware or by an improvement in the mathematical 
procedure. The hardware improvements, which are becoming available 
in the newer computers such as the IBM 360, are increased computa- 
tion speeds and faster access to external storage. 

Solution procedure improvement has been attempted by several 
methods. Some attempts are being made to take advantage of special 
structuring of the problem. These techniques are directed toward 
partitioning or decomposing the problems into manageable sub-problems. 
Another approach to improved solution time is the Primal-Dual 
algorithm which was developed by Dantzig. [3] 

This project has been directed towards another possible method 


for improving solution procedures. An appraisal is made of the step- 


by-step addition of constraints to the two-phase revised simplex pro- 
cedure as a possible solution technique which would serve to improve 
solution times. This method, called here the ''step-by-step addition 
of constraints (SSAC),"! exploits the advantage of obtaining a rapid 
solution to a small sub-problem and then moving from one optimal 
solution to an optimal solution of a slightly larger sub-problem as 
the remaining constraint equations are added one at a time. In the 
earlier stages of the solution procedure, the size of the matrices 
used in the iteration procedure would be relatively small, compared 
to the complete problem. By exploiting the advantage of multiplying 
smaller-dimensioned matrices, with the attendant shorter computa- 
tional times, a shorter overall solution time might be achieved in 
spite of the fact that a greater number of iterations would be 


required. 


2. Notation 

The notation used throughout this paper was chosen to correspond 
to the notation most frequently used in linear programming texts. 

(1) Upper case letters represent matrices. 

(2) Lower case letters represent columr or row vectors. 

(3) Subscripted lower case letters represent elements of row or 
column vectors. 

(4) Tableau notation, as illustrated below, is similar to that 


used by Dantzig. (3| 


z Ww Xe xX, Xe X3 Xb b 


o> 2G = 4 MR yt 





The columns corresponding to basis vectors are indicated by a dot. 
The pivot column (or row in the case of the dual simplex algorithm) 
1s indicated by an arrow. The pivot element, as determined by the 
appropriate minimum & criterion, is indicated by circling the element, 
and the basis variable which is then to be driven out is indicated by 
circling the associated dot. 

The tableau rows are referred to as the ''z'', "'w'' or Ri'' row where 
1 iS an integer corresponding to the sequence in which the constraints 
appear. That is, 'Rl'' refers to the first constraint, |"'R2"' refers to 
the second constraint and so on. 

Specific notation which will receive repeated use in this paper 
includes: 


m number of constraint equations 


number of variables 

mx n matrix of coefficients of the constraint equations, 
having elements aij 
mx m matrix of basis eons 

inverse of the basis 

the m-dimensional vectors which make up the B matrix 
an n-dimensional column vector having elements X 5 
the n-element row vector of cost coefficients having elements Cj 


the m-element column vector of the right-hand side of the con- 


straint equations (requirements vector), having elements b. 


3. Formulation of the Problem 


The general linear programming problem is stated as: 


Maximize 
ex 
subject to: 
Ax = b, 
and, 
x > 0. 


For the algorithm to be investigated, we define 
Z = "ee, 
We can then rewrite the problem as: 
Maximize 
2» 


subject to: 


and, 


To obtain an identity matrix to begin the two-phase revised 
simplex procedure we add artificial variables xXp415 Xpids oosey Xnim 
to the constraint equation with x,4j > 8 For gl "1 , ccm te 
If we define 
Wwe “X n+] “Xn4+2 meee Xnim 
and then maximize w, we will drive the artificial variables out of 
the basis and either obtain an initial basic feasible solution or an 
indication that the problem is infeasible. The process associated with 


maximizing wis usually referred to as Phase I. 


During the second phase (Phase II), the problem is stated in the 
following forms: 


Maximize 


subject to: 


411%] a 212% a @eee8e 1 ainkn I X n+] = by 
Pia y SUS 6500 ECD ne? my, 
am1*1 ate am? X2 + 96: 0s stem amn*n +X nm = Dm 
—C 1%] = C9XD —-— @eov000e ™ CH* A spe = 0 
Xt] tXnt2t se tXpem? W =a 


and, 
x; 20 Otte 15 250. veep. Me 
Because of the w equation we do not have an identity matrix to 
use aS a Starting basis. Therefore, we subtract from that equation 
each of the other equations, as appropriate, to remove the Xpj41,5 eee; 
Xntm Variables. We obtain the following form for the w equation: 
WH amt2, 1X1 * Ame2,2%2 t+ cee + Ame2 n%n = Pmr2 


where, 


m 
“mt2 y j = -s De enn i = 1% eeey NN 


and, 


m 
baad =O >a b. 


for those artificial variables, i, having non-zero prices. Our protlem 
iS now in a form such that we have an identity matrix to use as our 


Starting basis. 


Maximize 
a, 
subject tos 


sail ia late oe a4 949 + coe es «OF S1n*n eX 


| 9290) ne Pht? 2 
am1* q -+- amn2*o + seevseoce + Smn*n 5 Xn+m m 
—C 1%] 8 C4 X9 — e@eeoee ae Ch Xn Ts zZ 
me? , 1% 17 eme2 ,2%2 eet amt? -n*n + W 
and, 
eer for j = 162, SRL, nem 
If we let a - = =c. , the problem resolves into the three 
m1, j J 


matrices used for computations in the revised simplex procedure: 


a1] mB eoeeeeoeeeee se Fin 


a4 a99 @®eeoeoeeeeeees 29H 


p= 
an an? oeeoeoeeeees8e8 @ aan 
“itll: aan omtl yn 

Sit Duslumes:? * 2.0\0 «ime Omt2,n } 
| 0 oe@ece 8 @ 0 
] Os ° 

5 =" lk ; -}| , an identity matrix, and 

0 | 


L 


For simplification in handling the matrices when coded in FORTRAN, 
they will be re-arranged so that the objective function and the 
modified w equation appear in the first two rows rather than the last 
twos 


oeoeoeeeeeeeeoeee28e88 a 


ically m+1yn 


omt2, 1 eeeeveveveeeeeee @ amt2,n 


aa, @®eeeseeeoees21eoeesveeetee @ ay 


= A 
apy) eeoeeeeen7neeetet ee 8 8 oe SO 
am] @®eeeseeeseeeo2sv#e@eese#sesr?sesre*é ann 

0 
b+? 
ej by 
B = I, and b= : 
b 
m 


In the revised simplex procedure, the original A matrix and b vector 
of the full starting tableau (''original' tableau) are used at each 
iteration with the inverse of the current basis matrix to determine 


certain unknown elements of the current tableau. 
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When a constraint is added, the original A matrix and b vector are 
changed. The size of each is increased by one row, and the A matrix is 
increased by one column. The basis matrix is therefore increased in 
size by the addition of a row vector, § , corresponding to the elements 
of the new constraint which are in the current basis. It is, at the 
same time, increased by a column vector corresponding to a new artifi- 
cial variable. 

The new (m+1) x (m+1) basis is: 

B CG 
nt l 


By partitioning of the B matrix, the new inverse is readily 


B = 


obtained: 
ee 0 
~~ - ed 

Upon adding a new constraint to a problem, one of three cases will 
result: 

(1) the new constraint is satisfied; 

(2) the new constraint is not satisfied and the value of the 
additional artificial variable is positive in the basic solution; 

(3) the new constraint is not satisfied and the value of the 
additional artificial variable is negative in the basic solution. 

In the first case, the new constraint has no effect and the 
optimal solution to the entire original problem has been obtained if 
no further constraints are to be added. The optimal value of z is not 
affected by the new constraint. 

In the second case:, the two-phase revised simplex procedure is used 


to first drive out the artificial variable and then to maximize z. 
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In the third case, if the new artificial variable is assigned a zero 
cost in the z equation, the dual simplex procedure can then be applied to 
the infeasible primal to obtain an optimal solution to the dual, and 
hence an optimal solution to the primal. The artificial variable is con- 
sidered to be a legitimate variable of the original problem in the dual 
simplex approach. As a consequence, it never appears in the w equation. 
For convenience, we will carry along the 'w'' row of the tableau, for use 
when later constraints are added. 

When a solution is obtained which satisfies the added constraint, 
the optimal solution to the original problem has been found if there are 
no new constraints to be added, 

The addition of a new constraint to a linear programming will have 
the effect of either decreasing the previously obtained maximum solution 
or leaving it unchanged. That is, letting the subscript on z denote the 
number of constraints, 

Max Zoi, < max Z, . 

Solution of a linear programming problem by step-by-step addition 
of constraints may result in one or more unbounded solutions to the sub- 
problems if the initial constraints have fewer variables than are in- 
cluded in the objective function. This causes no difficulty, however, 
as the addition of one or more new constraints will serve to place 
bounds on the problem, unless the original problem is unbounded. 

The method of step-by-step addition of constraints has the advan- 
tage that an infeasible solution at an early stage will determine that 
the original problem has no feasible solution, and no new constraints 


need be added. The solution procedure is then terminated. 


4. Sample Problem 
Consider, aS an example, the problem: 


maximize 


subject to: 


and, 


x, 20 rer j = Wee, oe 


Adding slack variables and rewriting the problem, we have 


maximize 
23 
subject to: 
= 2x, = 4X 
3x, + 4x, a 6x4 4 


4x 7 3X9 


dn 
No 
~< 

Wd) 


Xp + Xp + 4x3 
and, 


Ig ee 39 


x; 20 for j 
For our solution)by the SSAC procedure, 
for the first sub-problem will be: 
W - 3x] - 4X - 6x3 - xX, = ~2u 


because the original w equation is 


wt x¢ = O, 


=) 
Aiy =e; 
ae 24 
=a 
5 


the modified w equation 


where X6 is the artificial variable introduced into ''R1''. 


The initial tableau will bes: 


Z W X6 X1 Xo X2 Xho Xe b 





In Phase I we will maximize w. We choose the most negative value in 
the modified w equation and pivot on the element which meets the. mini- 
mum ~@ criterion according to the usual (primal) simplex procedure. 
Note that in the first tableau the first pivot will always be such as 
to drive out the first artificial variable. 

We pivot using the product form of the inverse. [5| The | vector 


corresponding to the X3 column will be 


4/6 2/3 
N= 6/6 { =] 1 
1/6 1/6 
and, 
I Ome /3 
E = 0 l l 
0 o 1/6 
The new inverse is determined aa = EB as follows: 
I 0 2/3 I 0 0 I 0 2/3 
pe nn nn Cecilio ScmenaD 
0 oOo 1/6 & On 0 0 oOo 1/6 


The new tableau is now 





where X3 is now a basic variable; all values of Zy-C; > 0, and we have 
achieved the ''first!' optimal solution. 

To add a new constraint, the "'original'' tableau is augmented by 
one row and the column vector for X75 where X49 is the artificial variable 


associated with ''R2'', Ignoring the ''w'' row for the time being, the 


augmented "'original'' tableau will be 


z W X6 X7 XX, X2 X3 Xe XE b 





We compute the oy vector by determining from the previous solutions 
which vectors are in the basis. Observe that the coefficient vectors 
associated with z and w will always Be in the basis since we are maxi- 
mizing z and w. The first two elements of the & vector will therefore 
always be zeros. The third element in this case will be the coefficient 


of x. in the new constraint. 


3 
¥- | 0 12] 


The augmented inverse of the basis is determined by partitioning: 


The product - yao , in thapmease, is 


1 0 2/3 
0 0-12 | 1 = [o 0-2] , 
0 aye as 
and the new inverse becomes 
| Onen2/307 0 
wet 0 1 1 0 


Our next step is to determine if the value of the new artificial variable 
is positive or negative in the basic solution. We find 


X7 = 2u = 12x = ao. 


Because it is negative, our procedure tells us to assign c7 = 0 to the 
new artificial variable, X7, and use the dual simplex algorithm to drive 


out the artificial variable. 


1 9 x3 Xp X b 


0 87350 273...0 
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We pivot on the row having the most negative b. element in the dual 
simplex method. As a consequence, the new artificial variable will be 
dropped as a basic variable. 


Our pivot element is determined for ab j <0 


@ = Min [-2 15/245] = Min 0, 8/15, 2/6 | = 


The new inverse is determined by the product form of the inverse. 


0 | 0 0 0) 
0 0 | 0 0) 
| = — 
1/4 Om Oe 178 
1/2 by. oe. 


The current tableau becomes: 


Onis O.278 6 







0-1/3 1/4 


0 





it -Wal 
= Min | 32/7, 2| =2 

Since the up-dated requirements vector, b, still has a negative com- 
ponent, the dual simplex algorithm must again be applied. Pivoting on 


the xy, column gives: 


z Ww X¢ X 7 X1 Xo X3 Xu Xo b 





2 | 


Since all elements of b are now non-negative, our solution is optimal. 
Notice that the ''w'' row is retained but not operated on in the use of 
the dual simplex approach. 

Adding the third constraint, the augmented "'original'' tableau 


becomes 


XG aS) XQ xX] 





Again, the '"'w'' row is ignored until after we determine if the new 
artificial variable, XQ» iS positive or negative in the basic solution. 
After bringing in the new constraint, 
y = [0 0 0 | 

and, 

fr - |o 0 0-1/4 : 
The value of the new artificial variable, xg, for constraint '"'R3"', is 

SS tS er: 
Since the value of xg is positive in the basic solution, we can use 
Phase I of the revised simplex procedure to move to the optimal solution 
for the complete "original! problem. The augmented ''w'' equation in the 


"original'' tableau will be 


Z Ww X6 X7 XB Xp XQ XZ XW XE b 





a 


and the current tableau becomes 


Z Ww X6 x3 Xg Xy Xo X3 xy Xp b 


3 6 
(1) 2 
&= Min [6/3, 2/1] = 2 


Pivoting on the appropriate element of the column having the most nega- 





tive value in the ''w'' row, we bring X3 into the basis. The resulting 
tableau is then 


Zz Ww X6 X49 Xg X11 Xo X3 Xx Xp b 





Since all of the coefficients of the ''w'' row are not non-negative, 
we again pivot on the column having the most negative element, or in 
the case of ties, the left-handed one of the tied columns. Pivoting 
on the x, column, we bring x, into the basis. 

After this iteration, we find that all elements of the ''w'' row are 
zero and at the same time, all z; - Cj > 0 so we have completed both 


J 
Phase I and Phase II. The final tableau is: 


23 





optimal solution of the "original'’ problem: 


7 =a0 

ip oa 

X3 = 2 

Mine 12 
and all other amu. 


J 
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5. Programming Technique 

As a means of testing the feasibility of solving general linear 
programming problems by the SSAC method, the solution technique was 
coded in FORTRAN IVg for use on the IBM 360/67 computer. 

One subroutine was designed which would carry out the solution 
procedure by either the two-phase revised simplex method or the dual 
simplex method as appropriate. This same subroutine was used to solve 
the problems by both the standard revised simplex procedure and by the 
SSAC procedure. By using the same subroutine for both methods, it was 
hoped that any bias which might result from programming technique 
could be avoided. A driving routine was designed which would first 
solve a problem by the revised simplex method and then re-solve the 
same problem using the SSAC method. The two sections of the program 
were then timed. | Sections of the program not germane to the method 
being investigated, such as the reading in and printing out of data, 
were not included in the timing. The number of iterations required 


for solution by each method was tabulated. 


IThe timing routine was developed by Lt E. A. Singer, a student 
at the Naval Postgraduate School. 


ae 


6. Efficiency of the Algorithm \ 

several small problems, for which hand solutions could easily be 
obtained, were used for preliminary testing and debugging of the program. 
A number of larger problems were then solved to obtain a limited ex- 
perimental verification of the new procedure. The problems were chosen 
from three categories of problem types; mixing problems, transportation 
problems, (including a network problem and a transshipment problem), and 
caterer problems. Results based on this preliminary comparison, as 
Shown in Table I, were inconclusive. 

In order to obtain the solutions to a large number of problems, and 
aS a means of avoiding the considerable time and effort required to 
input data by hand, a routine was designed which would generate random 
problems. This routine utilized a random-number generator to generate 
elements for the A and b matrices. To insure the existence of a 
bounded optimal feasible solution, the problems were formulated as: 


Maximize 


subject to: 


and, 


x 2 0 ome = like. st 
a; D0 
b. >0 Ow 1) = 15 Seem. mM 


The problems were generated to have 70 variables, including slack 
variables, and 20 constraints. The distribution of the coefficients 


was uniform over the following intervals: 
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ay Uinintrerammos 1): 

b. , uniform (0,5); 

Cys uniform (-1,0). 

A total of forty-six problems were solved using this method of problem 
generation. A tabulation of solution times for these random problems 
iS given in Appendix I. 

A comparison of the solution times of these forty-six problems 
shows that the method of step-by-step addition of constraints was faster 
in thirty-four cases. The mean solution time for the SSAC procedure 
was 3.94 seconds faster than the mean solution time by the revised 
simplex procedure. By applying an appropriate statistical test to the 
solution results, it was determined that, with 95% confidence, the mean 
difference in solution times for the two methods is not less than’2.2/7 


seconds.” Therefore, it can be concluded that the method of SSAC is 


significantly faster than the revised simplex method. 


2 Detailed computations for the t-test and computation of lower 
confidence limits on the mean solution time difference are given in 
Appendix I. 


BAB LE! 


SOLUTION TIME RESULTS OF PROBLEMS USED IN PRELIMINARY INVESTIGATION 


Problem Revised Simplex Rddmewof Constr, Time Ratio 4 
(RS) (SSAC) | 
|} Iterations Time Iterations Time SSAC /E RES 
(sec. ) 


ixina Problems |} 






Waugh's a 9 . 106 a 
asoline yet 10 . 766 1.88 
ransportation Problems 

3 by 33 tableau 85 4118.829 043 

7 by 7 tableau S| 6.398 64 

3 by 5 tableau 18 6852 82 

3 by 4 tableau 13, | 446 Vaw7 

Teansahiatene™ 18 | 4.223 74 

Network F low 17 Breil a 

Caterer Problem 

Wardroom vaste 35 \Oavos 5) 

Hadley Napkei 10 ~745 bees 
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7. Concluding Remarks 

It 1S important to emphasize that the solution time is a function 
of the way the computer program is written. In the code employed, a 
full set of artificial variables was generated for each problem. An 
algorithm which takes advantage of existing slack variables for an 
initial feasible solution might well prove to be faster than the 
present program. 

Round-off error and exponent underflow can greatly affect the 
solution technique. The use of the product form of the inverse in the 
pivoting operation relieves this situation somewhat. In the experi- 
mental algorithm a routine was employed which set any element having 
an absolute value less than .0001 equal to zero. A better method might 
be to use double precision mode for computations and then allow values 
smaller than .0001 to be carried along in the solution. 

By using a single subroutine for the simplex iteration procedure 
in both methods of problem solution, it was hoped that any inconsistency 
due to programming technique could be kept to a minimum. That is, 
necessary computations for the iteration procedure were carried out in 
the same sequence for both methods of problem solution. 

It is recognized that the problems selected are not necessarily 
a representative sampling of linear programming problems. The manual 
input problems were selected primarily because they were large enough 
to allow the step-by-step addition of constraints to be demonstrated, 
yet small enough to handle conveniently as data inputs. The randomly 
generated problems were considered to be of a size which was large 
enough to effectively test the procedure, subject to available com- 


puter time. 


fa) 


The results of this limited number of tests indicate that the SSAC 
method is worthy of further investigation. This further effort could be 
directed in one of several areas. The present method should be applied 
to a large number of more diverse problems in order to obtain a better 
data base for verification of the results already obtained, and to deter- 
mine more accurately the advantage of this method over the revised 
simplex procedure. At the same time, it would be possible to determine 
some bounds of effectiveness of this procedure as to the size and struc- 
ture of problems. 

Modification of the step-by-step procedure might be attempted to 
take advantage of an existing basis in the original problem so as to 
require generation of fewer artificial variables. An attempt might 
also be made to modify the step-by-step addition of constraints pro- 


cedure for application to the primal-dual algorithm. 
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SEEENDIX I. 


Statistical Testing of Random Problem Solutions 
If it can be assumed that the solution times obtained by each 
method form two normal distributions, we can test to determine if the 
mean solution time by the addition of constraints method, PSSAC is 
less than the mean solution time by the revised simplex method, Prs. 
Letting X- be the solution time obtained by the revised simplex 
procedure, and Y; the solution time obtained by the addition of con- 
straints procedure for problem i, a pairwise comparison of the results 
of the two methods can be made. 
Consider the hypothesis that the mean solution time difference, 
B= Pres - PSSAc? 1S non-positive. That is, 
Hor =a C0, 
with the alternative hypothesis, 
Hj: pro. 
To show that the method of addition of constraints is faster, we must 
be able to reject Ho. 


If we format statistic for n = 46 samples, 





where, 


and, 


ee 
Ds Bs 


(| 


S - —_—c 


n 
ric | i= 


and then compare it against tabulated values of the cumulative t dis- 
tribution, we would reject the hypothesis if 
Z “eX »n-1 
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In this case we choose a 100(1 -X ) = 95% confidence level. The com- 


putation of the t statistic for n = 46, S = 6.74 gives, 


t = sete 7 Ge 397 
6.74 


The tabulated value for Yo 05) 46 =) 1. GCS Hee t > t we 
= 3 


Xa,n-1 
reject the hypothesis that the revised simplex method is faster. 

A lower confidence limit on the mean difference in solution times 
can be determined from the expression 


Prob }G- avr ¢ runes : 


Since with 95% confidence 


(D- y)%n ¢ ty 
S 


we can solve for the lower limit on Pp: 


Tewts: < : 
~~ 


For the given data, 


3.94 - (1.68)(6.74) < py 
(6.785) 


or, 
2.27 SAL : 
So we can say with 95% confidence that the mean difference in solution 
times is no less than 2.27 seconds. That is, the mean solution time 
obtained by the step-by-step addition of constraints is at least 2.27 
seconds less than the solution time obtained by the revised simplex 


procedure. 
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Prob. 
No. 





SIL I EVE I a RAEN i AS A RE RS ; : : = . 


LO CO ™~) Ov Fw hm — 


TABLE II 


SOLUTION TIME RESULTS OF RANDOMLY GENERATED PROBLEMS 


23.90 
18.19 
14.97 
28.25 
2.0 92a 
23530 
27.24 
12110 
13.47 
13.47 
536 IS 
14.50 
9.72 
13.48 
13.14 
10.76 
15.53 
19.62 
24.38 
19.16 
22.23 
21.93 
26.35 
23.65 
14.40 
9.72 
14.73 
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